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Abstract 

For real time severe weather and flash flood forecasting, accurate rainfleld nowcasting an hour ahead is invaluable. 
Numerical Weather Prediction models are not useful in this context, because they take about 6 hours to “run up” from a given 
set of initial conditions. An altemative is to use the random fields estimated by radars or satellites and employ the technique 
of stochastic forecasting to anticipate where the rainfrelds might go. A difficulty in using such images of rainfall is that the 
rainfrelds are contaminated not only by the naturally chaotic structure of turbulence, but also by the measurement noise from 
the sensor. To get a smoother field, which has the right amount of information to give good advection prediction capability 
without removing the relevant detail, is not a trivial matter and is the aim of the investigation presented here. 

A data-driven method for extracting information, at temporally predictable scales, from spatial rainfall data is described, 
which extends the Empirical Mode Decomposition (EMD) algorithm into two dimensions. The EMD technique is used here 
to separate spatial rainfall data into a sequence of high through to low frequency components. It has been suggested in the 
literature that the lower frequency components of spatial rainfall data exhibit greater temporal persistence than the higher 
frequency ones. This idea is explored here in the context of Empirical Mode Decomposition, to prepare rainfall data for 
nowcasts based on the temporal evolution of the lower frequency components. The presentation focuses on the 
implementation and development of the two-dimensional extension to the EMD algorithm and its application to radar rainfall 
data, as well as examining temporal persistence in the data at different spatial scales. 

Keywords: Rainfall Nowcasting, Empirical Mode Decomposition. 


1 Introduction 

Spatial rainfall data contain information at a broad range of spatial scales (Schertzer and Lovejoy, 1987; Harris et al., 2001; 
Pegram and Clothier, 2001). It has been suggested in the literature (Seed 2003; Tumer et al., 2004) that the lower frequency 
components exhibit more temporal persistence than the higher ones; this premise is used here to prepare the data for nowcasts 
based on the evolution of the lower frequency components of space-time rainfall sequences. Examination of the (radially 
averaged) power spectmm (Figure lb) derived from a typical instantaneous estimate of rainfall rate obtained by weather 
radar (Figure la) indicates that most of the power, hence potential for deterministic prediction in the context of nowcasting, is 
contained in the low frequency components. In this paper we focus on a data-driven technique to extract the high frequency 
(less persistent in time) modes as the first step in a rainfall nowcasting scheme. The technique employed is a two-dimensional 
(2-D space) generalization of the one-dimensional Empirical Mode Decomposition (EMD) technique introduced by Huang et 
al. (1998). In a single dimension, EMD analysis produces a set of Intrinsic Mode Functions (IMFs) that are very nearly 
orthogonal; in two dimensions a set of Intrinsic Mode Surfaces (IMSs) is produced with similar quasi-orthogonal properties. 
Two-dimensional EMD appears to have been first introduced by Linderhed (2002) in the context of image compression; the 
key contribution in this paper is to introduce the concept of 2-D EMD as a tool for the analysis of space-time rainfall data. 
This paper focuses on the implementation and development of the two-dimensional extension of the EMD algorithm in this 
context. 



Figure la. An observed convective rainfall field measured by S-Band weather radar at Bethlehem, South Africa 
(colour scale indicates instantaneous rain rate in mm/hr). The image is 100x200 with 1km 2 pixels. 
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Figure lb. Radially averaged power spectrum of instantaneous rainfall rate from typical radar rainfall 

data shown in Figure la. 

In the application presented here, the least persistent IMS (exhibiting the highest local spatial frequency and least amount of 
spatial correlation - hence nearly white noise) is computed and removed from the raw rainfall data leaving a residual 
composed of the more persistent low frequency structural components in the data. This process is equivalent to using a low- 
pass spatial filter, based on the observed properties of the data rather than the predefined basis fimctions used in traditional 
Fourier or Wavelet scale decompositions. In Sections 2 and 3, simple theoretical examples, showing the power of EMD in 
one and two dimensions, are presented as a “proof of concept’’ before applying the procedure to observed radar rainfall data 
from Bethlehem, South Africa (Section 4). These complement and extend the original presentation by Huang et al. (1998) 
and Flandrin et al. (2004). Computational aspects relating to image processing and surface fitting are covered in detail 
(Section 4) and conclusions are drawn in Section 5. 

2 Empirical Mode Decomposition in a single dimension 

The basic idea embodied in the EMD analysis, as introduced by Huang et al. (1998), is to allow for an adaptive and 
unsupervised representation of the intrinsic components of linear and non-linear signals based purely on the properties 
observed in the data without appealing to the concept of stationarity. As Huang et al. (1998) point out in their abstract: "This 
decomposition method is adaptive and therefore highly efficient. Since the decomposition is based on the local characteristic 
time scale of the data, it is applicable to nonlinear and non-stationary processes." 

Few sequences of observations of natural phenomena are long enough to test the hypothesis of stationarity and frequently, the 
phenomena are patently non-stationary. This tacitly applies in the measurement of rainfall at a point or in space-time because 
sequences of rain are interspersed with dry periods and during the raining periods, the variability of the intensity due to 
mixtures of rainfall type (stratiform, convective, frontal) confound the stationarity definition. The EMD algorithm copes with 
stationarity (or the lack of it) by ignoring the concept, embracing non-stationarity as a practical reality. For a fuller discussion 
of the genesis of these ideas, see the Introduction of Huang et al. (1998), who also heuristically demonstrate the implicit 
orthogonality of the sequences of Intrinsic Mode Functions (IMFs) defined by the EMD algorithm. 

In the application of the EMD algorithm, the possibly non-linear signal, which may exhibit varying amplitude and local 
frequency modulation, is linearly decomposed into a finite number of (zero mean) frequency and amplitude modulated 
signals, as well as a residual fiinction which exhibits a single extremum, is a monotonic trend or is simply a constant. 
Although EMD is a relatively new data analysis technique, its power and simplicity have encouraged its application in a 
myriad of fields. It is beyond the scope of this paper to give a complete review of the applications, however a few interesting 
examples are cited here to give the reader a feeling for the broad scope of applications. Chiew et al. (2005) examine the one- 
dimensional EMD of several annual streamflow time series to search for significant trends in the data, using bootstrapping to 
test the statistical significance of identified trends. The technique has been used extensively in the analysis of ocean wave 
data (Huang et al., 1999; Hwang et al., 2003) as well as in the analysis of polar ice cover (Gloersen and Huang, 2003). EMD 
has also been applied in the analysis of seismological data by Zhang et al. (2003) and has even been used to diagnose heart 
rate fiuctuations (Balocchi et al., 2004). 

2.1 Computing the one-dimensional EMD 

The EMD algorithm extracts the oscillatory mode that exhibits the highest local frequency from the data (“detail” in the 
Wavelet context or the result of a high-pass filter in Fourier analysis), leaving the remainder as a “residual” (“approximation” 
in Wavelet analysis). Successive applications of the algorithm on the sequence of residuals produce a complete 
decomposition of the data into a set of IMFs. The final residual is a constant, a monotone trend or a curve with a single 
extremum. 
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The EMD of a one-dimensional data set z(k) is obtained using the following procedure: 

1) Setr 0 (£) = z(k) and setz = 1. 

2) Identify all of the extrema (maxima and minima) in ( k ). 

3) Compute a maximal envelope, max ,_,(&) , by interpolating between the maxima found in step 2. Similarly compute the 
minimal envelope, min^, (k ). Cubic splines (as suggested by Huang et al., 1998) appear to be the most appropriate 
interpolation method for deriving these envelopes in one dimension (Flandrin et al., 2004). 

[" max^^fc) + mhq _,(&)] 

4) Compute the mean value function of the maximal and minimal envelopes (k) = - -- 


5) The estimate of the IMF is computed from IMF t (k) = r t _i(k) - m t _i(k) 

Each IMF is supposed to oscillate about a zero mean and in practice it is necessary to perform a “sifting” process by iterating 
steps 2-5 (setting r t _ x = IMF { before each iteration) until this is achieved. 

6) Once the IMF^ has a mean value that is sufficiently close to zero over the length of the data (defined by a stopping criterion 
within some predefined tolerance s) the residual r t (k) = r t _ x (k ) - IMFj(k) is computed. Alternatively the sifting procedure can 
be stopped when the difference in the standard deviation of successive estimates of/MFJ falls below a critical threshold 
(Huang et al., 1998). 

7) If the residual r t (k) is a constant or trend then stop; else increment i and return to step 2. 
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Figure 2. EMD based signal separation; all IMFs are plotted to the same vertical scale. Top panel is the combined 
signal; lower 3 panels are the decomposition which recaptures, almost exactly, the original components. 


Figure 2 shows the EMD of a composite data series (shown in the first panel) that is the summation of a sine wave, a 
triangular waveform and a slowly varying trend. The compact representation obtained by EMD extracts (almost perfectly - 
except near the ends) the three separate time series (shown in panels 2 to 4) that make up the composite signal, without 
resorting to Fourier or Wavelet techniques with restrictive assumptions about the form of the underlying oscillatory modes or 
basis functions. Figure 3 shows the analysis of the same data, using Wavelet decomposition. Here a fifth order Daubechies 
wavelet basis was (arbitrarily) chosen for illustration purposes; this choice of basis fimction may not be optimal for 
detrending but serves to demonstrate a typical decomposition. Seven levels of decomposition were required before the trend 
became apparent; this decomposition is clearly far less compact and physically meaningful than the EMD results in this case. 
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Figure 3: Wavelet based signal separation - The ‘data’ are the same as in Figure 2, the vertical scale has been 
compressed for a compact presentation. An arbitrarily chosen db5 wavelet basis has been used. 
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A similar decomposition analysis can be carried out using Fourier techniques. The Discrete Fourier approximation of a signal 
can be defined in terms of the Euler-Fourier coefficients (a 0 ,a k ,b k ^) with k = 1,2(Equation 1). The coefficients are all 

that are required to reconstruct the series and any signal can be well approximated (as long as it satisfies the Dirichlet 
conditions), provided m is sufficiently large. In equation 1, F(xj) is the Fourier approximation of the signal y- at each of the 

n discrete (evenly spaced) data points Xj . L is the range of values Xj over which the data set is assumed periodic. 


F(x;) = — + X {a k cos(27rkx j /L) + b k sin{2n;kx j /L)} (1) 

2 k= \ 

a k = 2 X cos(27 zkx; / L)y .• / n (2) 

j =1 

b k = 2 X sin(27rAx,- / L)y • /n (3) 

7=1 

Figures 4 and 5 show the result of decomposing the data using a finite Fourier series. Figure 4 shows the first 5 harmonics; 
while Figure 5 shows the series reconstruction by accumulating the lower harmonics up to m. Computing the Euler-Fourier 
coefficients provides a compact approximation of the original signal but fails to extract physically meaningful information. 
The ability to determine meaningful structural information is clearly important in a nowcasting context, which cannot be 
bound by the periodicity assumption implicit in Fourier methods. 
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Figure 4. Fourier based signal separation, the first 5 of 75 components - The dashed lines show the sine component 

and the solid lines the cosine component. 
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Figure 5. Reconstruction of the signal from the sine and cosine components, m represents the number of Euler- 

Fourier coefficients used in each reconstruction. 
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3 Empirical Mode Decomposition in two dimensions 

In two dimensions the EMD process is conceptually the same as for a single dimension, except that the curve fitting exercise 
becomes one of surface fitting and the identification of extrema becomes (a little) more complicated. Very little work appears 
to have been done which applies the EMD technique to two-dimensional data. Han et al. (2002) use EMDs in one dimension 
along four different directions to smooth Synthetic Aperture Radar (SAR) images and remove speckle. Nunes et al. (2003) 
develop a technique, which they term “Bidimensional Empirical Mode Decomposition” in the context of texture analysis in 
image data where they demonstrate several examples of intrinsic mode extraction from image data. Linderhed (2002; 2004) 
examined the use of EMD in two dimensions for image compression. Both of these implementations are very similar to what 
we propose in this paper. The 2-D EMD provides a truly two-dimensional analysis of the intrinsic oscillatory modes inherent 
in the data. Two-dimensional Fourier and Wavelet analyses are really applications of their one-dimensional counterparts in a 
number of principal directions. Fourier analysis concentrates on orthogonal “East-West” and “North-South” directions (e.g. 
Press et al., 1992). Wavelet analysis can, in general, consider any direction of the wavelet relative to the data, however a 
typical 2-D Wavelet analysis examines only horizontal, vertical and diagonal orthonormal wavelet basis fiinctions 
(Daubechies, 1992, pp 313; Kumar & Foufoula-Georgiou, 1993). In contrast, EMD produces a fully two-dimensional 
decomposition of the data, based purely on spatial relationships between the extrema, independent of the orientation of the 
coordinate system in which the data are viewed. 

3.1 Computing the two-dimensional EMD 

The algorithm follows intuitively from the one-dimensional case and may be briefly summarised as follows: 

1. Locate the extrema in the 2-D space including maximal and minimal plateaus. 

2. Generate the bounding envelopes using appropriate surface fitting techniques. We suggest conical Multiquadrics (for 
reasons explained in Section 3.2). 

3. Compute the mean surface function as the average value of the upper and lower envelopes. 

4. Determine the first estimate of an IMS by subtracting the mean surface from the data. 

5. Iterate until the IMS mean surface fiinction is close to zero everywhere. 

6. Estimate the IMS and Residual. 

7. If the Residual is a constant or a monotone trend, then stop; else return to step 2. 

3.2 Surface fitting for extremal envelope generation 

The generation of maximal and minimal envelopes is of key importance to a successful 2-D EMD implementation and is the 
most computationally intensive task. The problem is a familiar one of collocating a smooth surface to randomly scattered 
data points in two-dimensions. There are several options available to achieve this. Ultimately the fitting procedure reduces to 
computing the unknown value of the surface at a point s t = (Xj,yj) , by some linear (or nonlinear) weighting of the known 

data. In general, a basis function determines the influence of each known data point based on its spatial position relative to 
the unknown point s t . Nunes et al. (2003) use radial basis functions while Linderhed uses bi-cubic splines (Linderhed, 2002) 

and later chooses the more suitable option of Thin Plate Splines (Linderhed, 2004). We use radial basis fiinctions 
(technically, conical Multiquadrics), which are identical to Kriging (Cressie, 1991) with a purely linear semi-variogram 
model. It could perhaps be argued that it would be more appropriate to flt a semi-variogram model to the maxima and 
minima, but we feel this would be over-elaborate and presumptuous, as the extrema are only related by distance and cannot 
be considered drawn from a stationary correlated random field. Invoking Occam’s razor in the spirit of Huang’s original 
derivation of EMD, we wish to let the data do the talking and conical Multiquadrics assume the least structure of any linear 
surface fitting algorithm. 

The Ordinary Kriging estimate z t at any point i based on n observed data points is 


4 = I 4 z k 

k =1 



where z k are the observations and are weights associated with each observation and the target point. The mean is assumed 
unknown and the weights X k are constrained to sum to unity. The vector of weights X is obtained by solving the linear system 
in Equation 5 
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where y is a vector of semivariogram values, in this application simply defined by the linear distance basis fimction 


’u 


y(Sjj)= Sjj with Sjj the distance between point i and the y =1,2,...,« observation locations. T is the matrix of distances 


v 


between the observations, u is a vector of n ones and g is a Lagrange multiplier ensuring that the Kriging weights X k sum to 

unity, as required. The solution of Equation 5 is obtained using Singular Value Decomposition (SVD) in this application to 
ensure that a stable solution is assured (when the matrix is ill conditioned). This is achieved by truncating singular and near- 
singular components. 
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Although SVD is computationally less efficient than (for example) LU decomposition as a means of solving a dense linear 
system, it’s use is preferred here because of it’s robustness in the face of the near-singular Kriging systems which are 
frequently encountered in gridded data applications (Wesson & Pegram, 2004). 

A more efficient choice of interpolation technique would be useful and more work could be done in this regard, however care 
is required. Moving neighbourhood Kriging (a possible alternative to reduce the number of control points) can produce 
unwanted discontinuities in regions that are data sparse (Chiles & Delfiner, 1999, pp 201), such discontinuities would be 
amplified through the EMD sifting process. In addition, the particular choice of Ordinary Kriging as a method of generating 
the bounding envelopes was (partially) directed by the property that the estimates decay asymptotically to the mean of the 
observed extrema. 

3.3 Simpie two-dimensional EMDs 

In this section, applications of the 2-D EMD technique are presented. As an artificially constructed example Figure 6 shows 
the successful removal of noise added to a synthetically generated two-dimensional sine signal. The noise (with it’s high 
local spatial frequency) is almost completely described by the first IMS leaving a residual, which is closely representative of 
the underlying signal. 
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Figure 6. Example of EMD used for noise removal on a 2-D sine wave. The bulk of the additive white noise in the 

corrupted signal is well captured by the first IMS. 


Turning to a realistic example of the type we have been aiming for, Figure la showed an instantaneous radar rainfall field 
with an area of 100x200 km. A complete EMD of this field is shown in Figure 7 using a direct application of the 2-D EMD 
process described in Section 3.1; note the change in scale of the individual IMSs. The final residual (with a single extremum) 
gives a clear indication of the position of the largest convective raincell evident in the field. 
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Figure 7. Naive EMD of the observed rainfall field shown in Figure la - note the change in scale of the rain rates in 

the IMSs. 
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4 Application of 2-D Empirical Mode Decomposition to rainfall data 


The simple 2-D EMD application presented in the previous section is computationally burdensome when applied to rainfall 
data. In this section, to overcome this drawback, a number of specific refinements are presented which combine to make 
EMD tractable in practical real-time situations. 

4.1 Image processing techniques and optimisations 

Since an application of 2-D EMD requires the use of surface fitting techniques, large linear systems must be solved. The size 
of a system is determined by the number of known data points which are to be used in combination to find the unknown 
values of the surface at each remaining position in the field. The highly variable nature of rainfall data means that the field 
contains a large number of extrema from which the bounding envelopes must be constructed. Additionally there are a large 
number of zero (no rain) data, which constitute minima. By only considering raining areas, the size of the linear systems 
requiring a solution are greatly reduced since each raining area (if more than one exists) will contain a considerably smaller 
number of extrema than the entire data region and each can be treated separately. Furthermore, it makes no sense to consider 
an EMD in areas where the variable of interest does not exist, in this case the areas that are not raining. 

A number of well-known image processing techniques are implemented to isolate and process each raining area. Figure 8 
summarises the steps taken in processing the data with the boxes numbered 1-7 indicating different steps in the process. First 
a mask is generated to separate the raining and non-raining pixels (Figure 8, Box 1) in the instantaneous radar image; pixels 
below a threshold of 1 mm/hr are considered as non-raining and the remaining pixels are marked as raining. 





Figure 8: Summary of data processing; 1. Mask the wet and dry areas; 2. Trace the boundary of each wet region; 3. 
Separately label each wet region; 4. Decimate the fence by a factor of 5, isolating the ‘fence posts’; 5. Isolate the 
maxima in each sub-region; 6. Isolate the minima in each sub-region; 7. EMD analysis decomposes the data into the 
first IMS and the first residual using the maximal and minimal envelopes defined using the points in 4, 5 & 6. 


An outer boundary border-tracing algorithm (Sonka et al., 1999) is used to establish a boundary “fence” around each raining 
area (Figure 8, Box 2) and a flood-fill procedure is then used to fill each raining area with a unique identifier, resulting in 
separately labelled raining regions (Figure 8, Box 3). To reduce the computational burden of the algorithm even further, the 
boundary “fence” is decimated by a factor of 5 to reduce the continuous string of border points to “fence posts” while 
retaining the gross shape of the raining areas (Figure 8, Box 4). The next step in the processing of the data is to isolate the 
extrema in the rainfall field (Figure 8, Boxes 5 & 6). 

There are numerous possible techniques for identifying extrema in the rainfall fleld. Nunes et al. (2003) use a morphological 
reconstruction technique. One alternative, which was explored, is based on image segmentation and detection of extremal 
plateaus. However, our method of choice was to use a simple 8 neighbour search routine for identification of pixels with 
extreme values as done by Finderhed (2004). The choice was partly for convenience, but also because the majority of the 
(non-zero) extreme values in the rainfall fields studied tumed out to consist of single pixels. There is a rich literature on 
image processing techniques and the reader is referred to an introductory text such as Sonka et al. (1999) to explore the field 
further. Finally, the EMD analysis is carried out using the extrema within each raining area and the zeros at the “fence posts” 
of non-raining border pixels to specify the extremal envelopes (Figure 8, Box 7). Only one step of decomposition is shown 
here - the data is decomposed into the noisy first IMS and the first residual. 
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4.2 Results 

An analysis of over 800 individual radar scans, embodying mixtures of various ratios of Stratiform and Convective rainfall 
types, was carried out to determine the effectiveness of the 2-D EMD algorithm in separating the high frequency spatial 
components from the original rainfall data. Working on the basis that the average characteristics of the data over a range of 
spatial scales summarised by the power spectrum is intuitively useful, the (radially averaged) power spectra of (i) the original 
data, (ii) the first IMS and (iii) the first residual of each image were examined and compared. Figure 9 shows a typical result; 
the power spectrum of the residual shows a very close correspondence with that of the original data at high wavelengths 
while it contains far less power at the lower wavelengths. In contrast, the spectrum of the first (noisy) IMS has very little 
power relative to the data’s spectrum at high wavelengths but shows a strong correspondence at the lowest wavelengths. 
Figure 9 clearly indicates how the 2-D EMD technique moves the bulk of the high frequency components in the original data 
into the first IMS and leaves the high power, lower frequencies in the residual. Figure 10 shows a time average of this 
behaviour by plotting the mean values at each wavelength of the three spectra over five consecutive radar scans (beginning 
with the data used to produce figure 9). The radar scans are captured at approximately five-minute intervals. It is interesting 
to observe that the average of the spectra of the first IMS is flat for wavelengths longer than 10km, suggesting nearly white 
noise over this range. 
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Figure 9: Comparison of individual radially averaged power spectra of the radar rainfall data (of Figure la) with its 

EMD components: the llrst IMS and the first residual. 
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Figure 10: The same as Figure 9 but for the mean of individual power spectra for five consecutive, radar scans - 

Beginning with the spectra shown in Figure 9. 


The temporal persistence exhibited at the spatial scales represented in each of the three sequences of: (i) the data, (ii) the first 
IMS and (iii) the first residual was examined by considering their temporally consecutive power spectra. The notion of 
“spectral persistence” was used to determine how variable the spatial structure (at a particular spatial scale) is in time and 
hence to give an indication of the temporal predictive capability at each spatial scale. A summarised example of the analysis 
of a sequence of 5 radar rainfall images is presented in Figures 11,12 and 13 where a “matrix” of scatter plots is shown in 
each case. Scatter-plots of the pairs of power values at each discrete wavelength for five consecutive spectra (with the 1:1 
line indicated) are shown for, the original data (Figure 11) the first IMS (Figure 12) and the first residual (Figure 13). The 
rows and columns of the scatter-plot matrices are labelled from T 0 to T 4 and indicate separate radar scans between time T=0 
and time T=4. Each block in the scatter-plot matrix represents a scatter-plot of the power at each wavelength for the spectrum 
computed at T t versus that of the spectrum computed at T r 
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Clearly the plots on the “matrix” diagonal each compare a spectrum to itself and a perfect 1:1 relationship is observed in this 
case. For the off-diagonal plots, the degree of scatter amongst the data points indicates the degree of similarity between the 
spectra at individual wavelengths at increasing time lags with a large scatter indicating a weak similarity. The trends shown 
here are typical of the data analysed and show how the first (high average frequency) IMS has a temporally incoherent spatial 
structure, while the first (low average frequency) residual shows a temporally consistent structure. The behaviour shown in 
Figures 9-13 suggests that the high frequency IMS components in spatial rainfall data do not contain much predictive 
capability, supporting the suggestions of Seed (2003) and Tumer et al. (2004) to increase the degree of spatial smoothing and 
rely more on the information contained in the lower frequency components as forecast lead times increase. 



Figure 11. Spectral persistence scatter plots of the original data for a sequence of rainfall fields and those at 
successive intervals. This is constructed by plotting the values of power for each field at corresponding wavelengths 
coaxially. For example points A and B at the 10km wavelength are plotted against each other and appear ringed in the 
upper right diagram. 



Figure 12. Spectral persistence scatter plots of the sequence of l st IMSs of each pair of rainfall fields T 0 ,...,T 4 . 
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Figure 13: Spectral persistence scatter plots of the l st Residual for each pair of rainfall fields T 0 ,...,T 4 . 


5 Conclusion 

A new technique for analysing the spatial scaling structure of rainfall fields has been presented. The technique is a two 
dimensional extension of Empirical Mode Decomposition for the analysis of non-linear and non-stationary time series. An 
EMD analysis in two dimensions linearly decomposes the spatially distributed rainfall data into a set of Intrinsic Mode 
Surfaces (IMS), which are approximately mutually orthogonal and sum back to the original data. Each IMS contains an 
oscillatory mode inherent in the data at a different (narrow) range of spatial frequencies. The EMD analysis successively 
extracts the IMS with the highest local spatial frequencies in a recursive way, which is effectively a set of successive low- 
pass spatial filters based entirely on the properties exhibited by the data. The utility of the EMD technique for signal 
separation has been demonstrated in both one and two dimensions and applied to the analysis of a large set of 800 radar 
rainfall images in South Africa. The 2-D EMD technique is proposed here in the context of rainfall nowcasting to separate 
the less persistent high frequency components from the more persistent low frequency ones in the data. The aim is to remove 
the noisy high frequency components, which do not exhibit a strong temporal correlation and add little structural information 
to nowcasting algorithms. The scale separation achieved by 2-D EMD has been analysed using radially averaged power 
spectra to summarise the spatial structure of the data and filter outputs. In addition these power spectra have also been used to 
examine the temporal persistence of the spatial structure exhibited by the first IMS and it’s residual. The results presented in 
this paper agree with other work in the hydrometeorological literature, which suggests that the low frequency spatial 
components in rainfall data are most useful in a nowcasting context. This methodology is being exploited in ongoing research 
into rainfall nowcasting. 
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